Gust Alleviation using direct Gust Measurement 


Sven Marco Hoppe 

Flight Systems Research Center 
Department of Electrical Engineering 
University of California, Los Angeles 


March 10, 2000 


Contents 


1 Introduction 3 

2 Problem Definition 6 

2.1 System 8 

2.2 Observer 8 

3 Minimization of Accelerations 10 

3.1 Derivation of the optimal control law 12 

3.2 Necessary assumptions 17 

4 Minimization of Jerkiness 22 

4.1 Transformation of the problem to the previous case 24 

4.2 Necessary assumptions 26 

5 Computational Performance Requirements 28 

5.1 Minimization of accelerations 28 

5.2 Minimization of jerkiness 29 

5.3 Possible limitations 30 

6 Implementation with Matlab/Simulink 31 

6.1 Dryden model of gust 32 

6.2 Airplane state space model 32 

6.3 Controller 32 

6.3.1 Minimization of accelerations 33 

6.3.2 Minimization of jerkiness 35 

6.3.3 Bounds of the controller output 37 

6.3.4 Performance indices by means of H-infinity 40 

6.4 Evaluation module 43 

6.5 M-files 45 


1 


Contents 2 


7 Preliminary Simulation 46 

7.1 Test configuration, system and actuator 46 

7.2 Minimization of accelerations . . 49 

7.3 Minimization of jerkiness 53 

7.3.1 PT1 actuator 53 

7.3.2 PT2 actuator 58 

8 Application on a Wing-Fuselage System 63 

8.1 Test configuration, system and actuator 63 

8.2 Minimization of accelerations 70 

8.3 Minimization of jerkiness 76 

8.4 Discussion of appropriate sensor positions 83 

9 Final Remarks 85 

9.1 Conclusions 85 

9.2 Acknowledgments * • 86 

A Mat lab M-Files 88 

A.l const-def.m 88 

A. 2 con$t-calc-Lm 89 

A. 3 const-calc-2.m 91 

A. 4 Hinf-calc.m 92 

A. 5 pzmap^calc.m 94 


Chapter 1 

Introduction 


The increasing competition in the market of civil aircraft leads to operat- 
ing efficiency and passenger comfort being very important sales arguments. 
Continuous developments in jet propulsion technology helped to reduce en- 
ergy consumption, as well as noise and vibrations due to the engines. The 
main problem with respect to ride comfort is, however, the transmittance 
of accelerations and jerkiness imposed by atmospheric turbulence from the 
wings to the fuselage. This “gust” is also a design constraint: Light airplane 
structures help to save energy, but are more critical to resist the loads im- 
posed by turbulence. For both reasons, efficient gust alleviation is necessary 
to improve the performance of modern aircraft. 

Gust can be seen as a change in the angle of attack or as an additional 
varying vertical component of the headwind. The effect of gust can be very 
strong, since the same aerodynamic forces that keep the airplane flying are 
involved. Event though the frequency range of those changes is quite low, 
it is impossible for the pilot to alleviate gust manually. Besides, most of the 
time during the flight, the autopilot maintains course and the attitude of 
flight . Certainly, most autopilots should be capable of damping the roughest 
parts of turbulence, but they are unable to provide satisfactory results in 
that field. A promising extension should be the application of subsidiary 
control, where the inner (faster) control loop alleviates turbulence and the 
outer (slower) loop controls the attitude of flight. 

Besides the mentioned ride comfort, another reason for gust alleviation 
with respect to the fuselage is the sensibility of electrical devices to vibration 
and high values of acceleration. Many modern airplane designs - especially 
inherently instable military aircraft - are highly dependent on avionics. The 
fife time and the reliability of these systems is thus essential. To give a vivid 
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example, I would like to refer to a paper on Vibration Fatigue of Surface 
Mount Technology (SMT) Solder Joints [1] by S. Liguore and D. Followed. 
According to the graphs provided by the authors, the number of cycles to 
failure can be estimated by c ■ e~ p/ m , where p is the vibration level in Grms 
and c,m are positive constants. Figure 1.1 shows that if the vibrations are 
reduced by 50% the life time will increase by a factor 100-1000. 



Cycles to Failure 


Figure 1.1: Fatigue Life Correlation (Graph taken from [1]) 

Many papers propose controller designs based on structural measure- 
ment, e.g. by accelerometers, strain gauges or - more sophisticated - PZTs 
to detect gust by its effect. Some results show an improvement, but still, 
without more information, they will not be satisfactory. A very promising 
approach is the use of direct gust measurement. Then, information about 
turbulence is available before it has a significant effect on the airplane. As 
a result, the control system has more time for efficient countermeasures, 
especially if the sensors are located ahead of the wings. Investigation in 
that field already took place in the 1960s, e.g. for the F-100 Rough Rider 
turbulence measurement system [2]. 

Since then, the available measurement devices have undergone further 
development from simple mechanical vanes, as used for the F-100 exper- 
iments, to sophisticated laser/lidar systems. I would like to mention the 
research work of O.S. Alvarez-Salazar and G.M. Wang, A Novel Gust Moni- 
toring Device [3] which is based on the forward scattering of a laser beam, as 
well as the paper of D. Soreide, R.K. Bogue, L.J. Ehernberger and H. Bagley 
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on Coherent Lidar Turbulence Measurement for Gust Load Alleviation [4]. 

The reliability of these gust measurement data is essential for the perfor- 
mance of the alleviating system. Imprecise measurement can even make the 
controller worsen the effect of turbulence. Obviously, when the reading point 
is too far ahead of the wings, it is very likely that the turbulence changes 
between measurement and the encounter with the airplane. This problem 
is a matter of correlation time and will not be subject to my investigation. 
Besides, I believe that there is a practical borderline anyway, beyond which 
additional - even reliable - data have no further beneficial effect. 

Although the motivation for my work is the gust alleviation problem, 
I could also imagine further applications beside airplanes. In the field of 
automotive electronics, an active undercarriage for a sport utility vehicle is 
imaginable: Small radar sensors implemented in the front bumper detect 
the ground conditions and “road holes” in front of the tires and adapt the 
wheel suspensions for better ride comfort. 


Chapter 2 

Problem Definition 


In the course of the following, methods for controller design using data of 
direct gust measurement will be developed. Aside from a theoretical point 
of view the benefit of such measurement has to be weighed against the ad- 
ditional expenses incurred. 

In order to find the upper bound of possible improvements, a “perfect” sys- 
tem is considered. If the results are not satisfactory in this case, they won't 
be applied to a real airplane (with uncertainties about the model) at all. 

The following assumptions are made: 

L the whole system is exactly known, it can be described by linear dif- 
ferential equations 

2. the differential equations for structure and aerodynamics are at least 
second order 

3. the system is observable and controllable 

4. perfect accuracy of the measurement devices, particularly of the gust 
sensors 

5. no change of turbulence between measurement and encounter with the 
wings 

My investigation will cover both situations: 

1. sensor right at the location of the wing and its actuators (collocated 
system) 
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2. sensor ahead of the wings, e.g. abreast of the nose 

— > information about future gust is available ca. 0.1s before encounter 

The purpose of the controller is to reduce the effect of the gust encoun- 
tering the wing by providing suitable actuator signals for countermeasures. 
As opposed to many designs minimizing the displacements, I will deal with 
the accelerations, as the displacements (and also the velocity) themselves 
have no noticeable effect on the “interior” of the fuselage. Thus, I will con- 
sider the accelerations for the cost function as a first approach. Then I 
will also seize the suggestion of A. Tewari [5] who proposes optimal control 
considering the time-rate change of normal acceleration. He regards the 
“passenger /crew comfort or weapons aiming and delivery considerations” to 
be sensitive especially to this “jerkiness of the motion”. Both approaches 
seem applicable in view of electronic devices installed in the fuselage and 
the passenger comfort. 

As the differential equations of aerodynamic and structure for the gen- 
eralized displacements w are assumed to be at least second order, x can be 
defined as a state vector containing at least w and w. Therefore x includes 
the acceleration w. The time-rate change of normal acceleration (equals w) 
is not that easily accessible in any case (at least if I want to avoid z) but a 
sufficient approximation is possible. 

Hereafter I will make a discrete-time approach: The Riccati equation 
for a continuous- time control system is a serious problem, as there is no 
analytic solution for general system and weight matrices as well as general 
“gust functions”. In discrete-time a closed solution for a general case is 
possible. Necessary assumptions for the existence of this solution will be 
discussed. 

I will consider the location of the gust sensors variable to obtain a gener- 
ally valid result. Moving the measurement device to the nose of the airplane 
enables the engineer to get information about future gust encounters. Of 
course, modern lidar systems allow measurement far ahead of the airplane 
increasing the interval of prediction, however, when going too far, the re- 
liability of these data decreases due to jitter (uncertainty of time due to 
changes in aircraft speed etc.) and especially changes of turbulence before it 
encounters the airplane. For a theoretical analysis I consider the gust sensor 
data as absolutely reliable. I expect that the performance - starting from 
the measurement right at the wings - increases (the “cost” decrease) with 
the distance from the wings very fast at the beginning, but then will come 
to a saturation, since the knowledge of gust far in the future will certainly 
not have a remarkable beneficial effect on the control problem. I would like 
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to point out that - within this investigation - the expected behavior will not 
be due to the increasing uncertainty (as discussed above) of future gust data 
which I assume to be certain. However, for a real application the reliability 
of measurement is a constraint for the sensor position. The interval between 
measurement and encounter is r, so that z(i) is known for t < t < t + t. 


2.1 System 

The system can be described as a state space model in continuous-time. The 

matrices are assumed to be known exactly. 

x = Ax + Bu + Gz (2.1) 

y = Cx + Du (2.2) 

z(f) gust (velocity of vertical wind) directly at the wing. General case: 
Gust is represented by several reading points (where also the sensors 
are located) and a truly known shape function included in G. For the 
easier case of (in spanwise direction) homogeneous gust, z is simply a 
scalar. 

y(t) measurement vector of structural displacements (plunge, slope, torsion 
etc.) 

u (0 controller output 


2.2 Observer 

In view of the assumption that the model and the process noise (gust) are 
perfectly known for current time as well as the fact that there is no measure- 
ment noise, a Kalman filter is not required. A standard Luenberger observer 
is applicable: 

x = Ax + Bu + Gz + K (y — Cx) 

If (A - KC) leads to a stable system with K sufficiently high, the initial 
error will approximate zero after a short time with the limit: 

lim (x — x) — 0 

t-+o c v 

The observer will thus lead to almost perfect estimations soon after the 
beginning of the experiment. As in this case, a computer simulation is 
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processed, the state vector is accessible anyway. To keep the equations more 
simple and to focus on the controller design itself I leave out the observer 
(upper boundary of performance). 


Chapter 3 

Minimization of 
Accelerations 


In view of the underlying differential equations, I assume that the system 
is provided for continuous-time in state space, that is the matrices A, B, 
G are given. The aim is to minimize the time rate of acceleration of the 
systems, especially of the fuselage. I suppose that the sampling rate is suf- 
ficiently high, so that my goal is achieved by minimizing the acceleration at 
the sample times. 

The current time is /*., thus the purpose of the calculation is to find an opti- 
mal control output u k To insure that all available information is considered 
for the control output, the whole calculation leading to u k will be made for 
every (discrete) time t k . The rest of the sequence of control outputs will 
not be used. State of the art signal processors should be able to provide the 
result without a remarkable delay. The discrete-time approach easily allows 
to keep the sensor position variable. 

An optimal control approach for discrete-time systems is described in de- 
tail by K. Ogata [6] for a system without a (known) disturbance and using 
the state vector itself for the cost function. Due to these two differences and 
the fact that there are various cross terms within the cost function, I have 
to change the transformations resulting in another kind of Riccati equation. 


Cost function: 

J 

Xk 


2 X k+ N Sxk+ N + 2 [x^QiXi + u^QaUi 

i— k 

x(/)| (=ffc same notation for and z k 


( 3 . 1 ) 
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x k = *(<)U = (Ax(t) + Bu(«) + Gz(f))| t=tfc 

= Ax k + Bilk + Gzk 

x k +i = *x k + Tu k 4- ©z k 

* = e AT 

T = J\ At <1tB = A" 1 [e AT -l]B 

© = J T e Ar dTG = A " 1 [e Ar -l] G 
T = t k +i - tk V/c > 0 sampling period 


(3.2) 

(3-3) 


JV, or to be precise AT, is the length of the considered finite horizon. 


As mentioned above, I can assume that x contains w and w where w 
represents the vector of velocities. With regard to my approach I do not 
want to consider velocities within my cost functions and would choose Qi 
in such a way that only the terms dealing with w are unequal to zero. 
However, as for the following calculations, I will not go on such restrictions, 
especially not when discussing the necessary assumptions. 

Depending on the location of the gust sensors and on the horizontal 
velocity of the airplane, I will have (for a fixed T) a different number n of 
gust data z k . . . z k+n -i available referring to the two cases as described in 
chapter 2: 

1. gust sensor right at the location of the wing and its actuators, that is 
n = 1, no future gust data 

2. gust sensor ahead of the wings. The larger n > 1 is, the more is the 
measurement device moved towards the nose. I assume that n < N 

For the case of n > 1, it is necessary to know the time delay r between 
measurement and effect of the turbulence. Therefore, I propose to use a 
sensor right at the wings in addition. Thus, r - and finally /? - can be 
calculated as the maximum of the correlation between the sensor signals. 
Each measured value can be written in a shift register of length N where 
the variable wTite position is n. 

Of course, the question arises what values for the unknowm gust z k+n . . . 
z k +N-i will be used within the cost function. Obviously, gust can be seen 
as a process noise (unbiased filtered w^hite noise) for these times. For this 
additive-disturbance stochastic problem, the optimal control is identical to 
the deterministic case, that is z\ = 0 for i = k + n . . . A; + N — 1 (certainty- 
equivalence principle) [7]. The reason why I do not just set N = n is that 
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even if I don’t have information about gust for these times, I have to consider 
the dynamic of the system to suppress accelerations (e.g. due to oscillations). 
I will derive the optimal control law for somehow known gust z k . . . z k+ N_i 
and will later, when it comes to applying the controller, set the unknown 
parts to zero as described above. This method is possible since there are no 
partial derivatives with respect to Zj involved hereafter. In doing so I can 
avoid unpractical fall differentiations. The current state vector x k is known. 

3.1 Derivation of the optimal control law 

Using the Lagrange’s method, equation (3.3) provides the necessary con- 
straints: 

j jfc+JV-l 

L = -x£ +N Sx k+N + - Y [*iQi*i + 

** i=k 

+ ui r Q 2 u i + xf +1 (-X 1+1 -I- $Xi + Tuj + ©Zi)] (3.4) 

1 T o 

L = -x k+N Sx k+N + 

+ 7 : Y (Ax i + Bu i + Gz i ) r Qi(Ax i + Bu i + Gz i ) + 

* » L 
i—k 

+ uj r Q 2 Ui + A^ (-x i+ i + $x; + Tui + ©zi)] (3.5) 

1 T „ 

L = -x k+N Sx k+N + 

fc+JV-l 

+ - Y, [x^A^^QiAxi + x?A T QiBui + x- r A T QiGzi + 

2 i=k 

+ u?B T QiAxi + u^B T QiBui + u?B T QiGzi + 

+ z^G^^QiAxi + z^G t QiBui + z^G T QiGzi + u?Q2Uj + 

+ A(+i (— x;+i 4- * Xi + Tui + ©z;)j 
, , k+N-l 

L = 2 X k+N Sx k+N+2 £ [( U T Z T)- 

( A t QiA A t QiB A t QiG \ /x|\ 

B t QxA B t Q x B + Q 2 B t QiG Uj + 

\ G t QiA G t QiB G t QiG / \ ^ / 

+ A-^! (-x i+ i + ^Xj + Tuj + ©Zi)] 
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Using these substitutions: 


2 X k+N^ x k+N + 

, Jfe+N-1 . 

+5 £ (*r «. T ■?) 

■ Sn 

s 12 

Sis 

S 2 1 

s 22 

S 23 

i=k 

S 31 

s 32 

s 33 

+ Xj +1 (-X|+i + <&Xi -1- Tuj + ©zi)| 





Note: Qi = Qj and Q 2 = Qj > thus S? = Sji Vi, j S T = S 
Q 1? Q 2 and S are real symmetric weight matrices. 


The necessary conditions for L to have a minimum read: 

— — == 0 : x? Sn + uf S 2 X + zj S 31 + — A^ = 0 (3.6) 

OXi 

for i — k + 1 . . . k + N — 1 

x k-fN^ - ^k+N — ® 

for i — k + N 

|^ = 0 : x? , S 12 + u?'S 22 + z?'S 3 2 + A?; i r = 0 (3.8) 

CJlli 

for i = k . . . k + N — 1 

• 77 ^ = 0 : -Xj + ^Xi_ x + ruj_i + ©Zi_i = 0 (3.9) 

0 Ai 

for i = k + 1 . . . k + N 
Transformations that will be used later: 


Ai = 4> T Aj +1 + S 13 Zi + Si 2 Ui + SnXi 
for i = k + l...k + N — 1 
x i+ i = $Xj + Tuj + ©Zj 


(3.10) 

(3.11) 


The set of equations represents a two-point boundary-value problem: x k is 
known for first boundary and (3.7) provides information about the last. In 
order to solve this problem a kind of Riccati transformation is applicable. I 
would like to point out again that the Riccati equation for the discrete-time 
case is much easier to handle than the one for continuous- time and a closed 
solution is possible. Assuming that Aj can be written as follows: 

where "Pj = Pi 


Aj — P jXj + p; 


(3.12) 
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Taking (3.6) as my starting point: 

0 = x^Sn + u^S2i + z^Ssi + 

o = x?’s 11 ^- 1 r + u?'S 2 i#- 1 r + z?'S 3 i^- 1 r + 

+ a^t 

-xJ r S 1 2-u' r Ss2-zj r S32 

A| r ^ -1 r = X? (Su*-^ - S 12 ) + XI? (s 21 *- x r - s 22 ) + 

+ zj - S*a) 

for i = A; + 1 . . . k + N — 1 

and thus: 

uj = A^- 1 r(s 21 ^- 1 r-S 22) _1 - 

- X T (Sn$ _1 r - Sia) (s 2 i^ _1 r - S 22) -1 - 

- z 7 (S 31 *- l r - s 32 ) (Sai*-^ - S 22 ) _1 
Ui = (r T $- T s 12 - s 22 ) _1 r T ^- T Ai - 

- (r T *~ T S 12 - S 22 ) -1 (r T *- T Sn - Sai) x; - 

- (r T ^- T Si 2 — s 22 ^ (r T *- T s 13 — S 23 ) Zi 

' V ' 

Mjnv 

Ui = M inv r T ^- T A?' - M inv (r T ^- T Sn - S 21 ) Xi - 

- M inv (r T $- T Si 3 - S 23 ) Zi (3.13) 

for i = k + l...k + N— 1 

Substitution of (3.12) into (3.10): 

PiXi + pi = 4^ T Pi + i Xj + i +4> T pi+i + Si 3 Zi + S 12 Uj + SijXi 

4>xi+ru[+©zi 

PjXi + Pi = ^ T P i+ i^ + Snj Xi + ^ T P i+ i© + Sis^ Zj + 

+ 3> T Pi+l + (* T P i+1 r + S 12 ) Ui (3.14) 

for i = A’ + 1 . . . k + N — 1 

Substitution of (3.13) into (3.14): 

PiXi + Pi= # T Pi + i$ + Su- 
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- (* T P i+ ir + s 12 ) M inv (r T *- T Su - s 2 i)] xi + 

+ [* t P 1+1 © + S x3 - 

- (^ T P i+ ir + Si 2 ) Mj nv (r T ^- T s 13 - s 23 )] Zi + 

+ ^ T Pi+1 + ($ T p i+1 r + s 12 ) M inv r T #- T (P,x, + Pi) 

[-Pi + * T P i+ i* + s u - (# T p i+ ir + s 12 ) M inv • 

• ((r T *- T s n - s 21 ) - r T ^- T Pi)] Xi + 

+ [^ T p i+ i© + s x3 - (# T P i+ ir + Si 2 ) M inv • 

• (r T ^" T Si 3 - Saa)] z i + ^ T Pi+i + 

+ [(* T P i+ ir + s X2 ) M inv r T ^- T - 1] Pi = 0 (3. 15) 

for i = k + 1 . . . k + N — 1 

Equation (3. 15) must be satisfied for all xi, i = k + 1 . . . k + N - 1 
=> set term in first squared brackets to zero: 

-Pi + 4> T P i+1 4> + s n + * T p i+1 rM inv r T *- T Pi - 

- * T p i+1 rM inv (r T *- T Sn - s 2x ) + 

+ s 12 M inv r T *- T Pi - s 12 M inv (r T ^- T Sn - s 21 ) - o 
[1 - (* T p i+1 r + s 12 ) M inv r T $- T ] Pi = 

* T P i+ x (* - rM inv (r T $- T S u - S 2 i)) + Sn - 

- S 12 M inv (r T 4>- T Sn - Sai) 

Pi = [1 - (^ T Pi+ir + Sia) M inv r T *- T ] _I • 

■ '* T P i+ i (* - rM inv (r T ^- T Sn - S 21 )) + Sn - 

- S 12 M inv (r T ^- T Sn - Sax)] (3.16) 

for i — k + 1 . . * A: + N — 1 

Equation (3.16) is a kind of Riccati equation, that allows a backward calcu- 
lation of Pi for i = k + N — 1 . . . A: + 1. The necessary final value P k +N can 
be obtained by substituting (3.12) into (3.7): 

Pk+N x k+N + Pk+N = SXk+N 
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Since there is no explicit constraint for x k + n (this is not applicable in view 
of a permanent disturbance by the gust), this equation must be satisfied for 
all x k+ N. As a result: 

Pk+N = S (3-17) 

Pk+N = 0 (3.18) 

Now, Pjc+N • ■ • Pk+i are known and can be used to calculate Pk+N-i • • - Pk+i 
For this reason, I would like to remind of (3.15). The last two rows must 
equal zero as well: 

[* T Pi+i© + S i3 - (* T P i+ ir + Sia) M inv (r T *- T Si 3 - S 23 )] *i + 
+ * T Pi+1 + [(^ T P i+1 r + Si 2 ) M inv r T *- T - 1] Pi = o 

Solving for pi: 

Pi = [i - (^Pi+ir + s 12 ) M inv r T *- T ] -1 • 

• [* T Pi+1 + (<I> T P i+1 © + s 13 - 
- (* T P i+ ir + s 12 ) Mi nv (r T *- T s 13 - s 23 ))zi] ( 3 . 19 ) 

for i = k + 1 . . . k + N - 1 

As Zk • • Zk+N-i are (somehow) known, pk+N-i • ■ • Pk+i are available. 

When it comes to calculating the required control output Uk, equation (3.13) 
is not applicable, since it only holds for i = k + l...k + N — 1. Therefore I 
refer to (3.8): 

0 = x^ S 12 4- uj S 22 + Z\ ^32 + iT 

0 = X? S12 + uf S22 + Z T ^32 + (xT +1 p i +1 + pf+l) r 

0 = xf S12 + uj S22 + zf S32 4 - ((*Xi 4 - Tuj + ©z 0 T Pi + i 4 - Pi+i) r 

0 = x^ ^Sx2 + $ T P i+ 1 r) + u? ^22 4- r T p i+ ir) 4- 

+ zj (S32 + © T Pi-hir) + p{^xT 

0 = (s 2 x + r T p i+ x^)x i + (s 22 + r T p i+ 1 r)ui + 

+ (s 23 + r T p i+1 ©) Zi + r T P i+! (3.20) 

for ?' = A: ... k + N — 1 

uj = — (S 22 4 - r T Pi_j-xrj ^S2x + r T Pi + i^j xi + 

4- ^23 + r T p i+1 ©) zi 4- r T Pi+1 
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Finally the optimal control law reads: 

u k = -(s 22 + r T p k+ 1 r ) _1 [(s 21 + r T P k+ i*)x k + 

+ (s 2 3 + r T p k+1 ©) z k + r T p k +i] (3.2i) 

To compare the performance of the gust alleviation system a “random” 
gust function, or to be precise a sequence of data, should be previously 
recorded for a period of time to guarantee equal conditions for the tests. 
Then, for different values of 7) and also of Qi, Q 2 and S a given system 
(represented by the matrices A, B, G) will be controlled. In order to obtain 
comparable results I introduce a set of cost functions: 

J x = f x T (f)RCx(f)df ride comfort (3.22) 

J to 

J 2 — [ x T (f)AD x.(t)dt “aggregated” displacements (3.23) 

Jto 

J 3 = [ u T (t)CE u (t.)dt control energy (3-24) 

Jto 

The state vectors used for these functions refer to the case of the controller 
applied on the continuous-time model as described in (2.1). This is closer to 
a real application, in which case the accelerations between t-i and t,+i have 
to be taken into account as well. I would like to point out that Qi , Q 2 and S 
are not applicable to J\, J 2 and J 3 as they only establish priorities within the 
controller design and may very. For comparable results, the weight matrices 
need to be constant throughout the tests. 

3.2 Necessary assumptions 

The above-mentioned approach is only applicable if and only if all of the 
following matrices are nonsingular: 

A 

A -1 exists if and only if det(A) ^ 0. In view of an observable and 
controllable system this should be the case. At least if I assume that 
the system was already stabilized by a slower feedback controller (the 
main purpose of the controller designed above is the minimization of 
accelerations) the eigenvalues of A can then be assumed to be all neg- 
ative. As a result, A would be positive definite and thus nonsingular. 
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S 2 2 


S 22 — B t QiB + Q 2 

If (B t QiB + Q 2 ) is a positive definite matrix then its determinant is 
unequal to zero. 

If a matrix G is posit ive definite, that is J = v T Gv > 0 V v and J — 0 
only for the case of v = 0, then J = (Ty) T G(Ty) = y T (T T GT)y = 
y T Gy is a positive definite form, too. 

If A is positive definite and B is at least positive semidefinite, then 
v t Av -f v t Bv = v t (A + B)v is a positive definite form, so that 
(A + B) is a positive definite matrix. 

If the terms of Qi that refer to nonzero parts of B are at least positive 
semidefinite then B t QiB will be, too. Assuming that Q 2 is positive 
definite, exists. 


$ 




„A T 


OO 


= £ 


(AT)* 

A;! 


It can be shown [6] that the sum converges absolutely for all finite T, 
so that exists. Therefore, the discrete time approach is possible for 
any A 

According to [6], e AT is nonsingular and 4> -1 = e~ AT . 


(s 2 i*- 1 r - s 22 ) 

I have to prove that the determinant of (S 2 i‘I> -1 r — S 22 ) is unequal 
to zero. I will use | | equivalent to det(). 


|s 21 *- 1 r - s 22 | 

B t Qi A# _1 A _1 [^-I]B-B t QiB-Q 2 | 


/0 

io 


Interim calculation: 

A^ -1 A _1 = A 


I+A( _ r) + A!tl£ + A!fcl)! + . 


2! 


3! 
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A + A 2 (-T) + 


a 3 (-r) 2 a 4 (-t) 


2! 


+ 


3! 


+ 


1 + AhT) + ^z£l + ^zDl + . 


2! 


3! 


= 

The simplified equation reads: 

B t Qi^ -1 [# — I] B - B t QiB - Q 2 | ^0 
B t QiB-B t Qi^- 1 B-B t QiB-Q 2 | # 0 
!b t Qi^- 1 B + Q 2 | 7^0 


-1 


If (B t Qi 3> -1 B + Q 2 ) is a positive definite matrix its determinant is 
unequal to zero. This is the case if one addend is positive definite and 
the other one is at least positive semidefinite. Q 2 is assumed to be 
positive definite (see S 22 ). 

If (Qi$ -1 ) is a positive definite matrix, (B t Qi*I> - 1 B) will be, too, 
so that the sum is positive definite. However, I think the criteria of 
positive definite property is too restrictive in this case. 

Simplifying the expression: 


B t Qi$ _1 B + Q 2 
IQ2||i + Q^ 1 b t Qi^ _1 b 


7^0 

#0 


IQal + 0 as Q 2 is positive definite, so that the resulting condition 
reads: 


I + Q2 1 B t Q 1 #- 1 B 


# 0 


Using the following relationship: 

|I„ + AB| 

A ~ (n,r)-matrix, 


|I r + BA| 

B = (r,n)- matrix 
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Then: ( 

|l + # _1 BQ2 1 B t Qi| 7^ 0 

Because of det(^) 7 ^ 0 I can multiply with ^ 


^ + BQ2 1 B t Q 1 


7^0 


As Qi and Q2 are almost free parameters, it should be possible to 
fulfil this condition, so that (S 2 i* -1 r - S 22 ) _1 = M? v (and thus 
also Mj„ v ) exists. 


[1 - (* T p i+1 r + s 12 ) M inv r T *- T ] 

Again, I have to prove that the determinant is unequal to zero. Using 
the relationship |I n 4- AB| = |I r + BA|: 

i - (* T p i+1 r + s 12 ) M inv r T *- T | ? o 
i - M inv r T *- T (* T p i+1 r + s 12 ) | ? o 

i - (r T *- T s 12 - s 22 ) _1 r T *- T (* T p i+1 r + s 12 ) | # o 
(r T *- T s 12 -s 22 ) _1 |- 

V ' 

• r T ^- T Si 2 — s 22 — r T $~ T (^ T p i+ ir + Si 2 ^ | ^ o 


s 22 + r T #- T ^ T P i+ ir 


s 22 + r T p i+1 r 


7^0 


S 22 was chosen to be positive definite above. For completely state 
controllable systems, it can be shown that Pj ( i is positive definite 
or positive semidefinite [ 6 ]. The sum of a positive definite and an 
at least positive semidefinite matrix is positive definite. Therefore 
S 22 + r T P i+ ir is unequal to zero and the matrix 

I - ($ T P i+ ir + S 12 ) M hlv r T *- T ] exists. 
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(s 22 + r T p k+ 1 r) 

Nonsingular as S 22 is positive definite and P k +i is at least positive 
semidefinite (see above). 



Chapter 4 


Minimization of Jerkiness 


The idea of this approach is that time rate changes of acceleration - the 
jerkiness - have a direct bearing on the passengers comfort and on the 
strains for electronic devices mounted in the fuselage which are delicate to 
such vibrations. In simple words: Fast changes of acceleration are considered 
to be worse than even larger but almost constant accelerations. However, 
as already mentioned above, the time-rate change of acceleration is not easy 
to handle. To illustrate this problem, I would like to give an example: 

y = w 

w = a\w + a 2 W + z + u 

y — aiw + ayw + z + ii 

If I want to minimize the jerk V of the system using u provided by a digital 
controller I will run into trouble since u(t) is a “staircase function’ 1 which 
is in principle not differentiable for all t = tk and will lead to high peaks 
for these points in time. From a pure mathematical point of view I could 
use the generalized derivative that deals with 6-functions for these times, 
but one the one hand I would obtain an integral of quadratic terms of 6(t) 
within my cost function which is really not easy to handle and on the other 
hand it would remain an ill-posed problem since any change of the controller 
output leads to new and even higher values of jerkiness than the gust (as a 
kind of lowpass filtered noise) itself. 

More generally speaking, the difference of order - the first derivative of y 
with respect to time on which a has an instantaneous influence - between 
input u and output y of the whole system needs to be at least three. If, for 
example, the differential equations of structure and aerodynamics are second 
order, the used actuator or a smoothing filter for the controller output must 
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be at least of first order. 

The structure of A and x is not prescribed. To maintain this generality 
I will use the following substitutions: 



where w2 is the vector of generalized displacements for which I want to min- 
imize the jerk. Note: wl refers to the actorator(s), see the m-file consLdef.m 
at the appendix (A.l). With respect to the pre-condition of order, u has no 
instantaneous effect on w2. Thus 


w2 = Rx + G2z 


The cost function for the approach of minimizing the time-rate change of 
acceleration reads: 


fc+JV-l 


1 Y K *r iV " 

J — Q X k+r>jSxk+N + 9 52 


i—k 


vi2j Qiw2j + Ui Q 2 uj 


(4.1) 


Note: \v2j = w2(t)| f 

Inserting the substitutions into the cost function: 


k+N-l 


1 1 * 1 r 

J = 2 X k+N^ x k+N+ 2 ® Zi) T Qi (Rxi -f G2 Zi) + uj Q2Uj 


i=k 

The corresponding Lagrange’s function reads: 

fe+iV-l, 

2 " ‘K+l’N ivti^ ' 2 


1 1 

L = + 2 5Z ^(Rxi + G2 zj)^ Qi (Rxj + G2 i\) + 

i~k 


+ Ui r Q 2 u i + (-X|+i + ^Xi -I- Tui 4- ©Zi)] 


(4.2) 
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4.1 Transformation of the problem to the previous 
case 

It would be a convenience to benefit from the calculations and finally from 
the control law derived in the previous chapter. This is in fact possible as I 
will show hereafter. 

First of all, I have to find an approximation for i\ since this information is 
not available. Knowing that gust can be seen as lowpass filtered white noise, 
I can assume that z (t) is smooth and, especially, does not jump. Therefore 
the approximation 

Zi « j; 

with T sufficiently small is applicable. 

Then I will introduce 


Zi = 


Zi 

Zi-l 


so that I can express i\ in terms of i\\ 

Zi = i[l -I ] Si 

To avoid dealing with two variables for the gust I will also rewrite the system 
equations: 


x = Ax + Bu+|^G ojz 
x k+1 = $x k + ru k + [ © o J z k 


(4.3) 

(4.4) 


Using these substitutions within the Lagrange’s cost function (4.2): 


L = -x£ +N Sx k+N + 

fc+iV— 1 


+ 5 E 


RAxi + RBui + R [ G O ] z s + G2 ^ [ I -I ] z, 
■ Qi (rA Xi + RBui + R [ G O ] z- t + G2 -t [ I -I ] Zj) + 

+ Ui r Q 2 u i + (-x i+ i + ^Xi + Tui + ©zj) 
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The last substitution I will make combines the coefficients for zj: 

G = R [ G O ] + G2 ^ [ I -i] 

This simplifies the Lagrange’s function: 

1 1 k+N— 1 r ~ \T 

L = -x£ +N Sx k+N + - S (RAxi + RBui + Gi i) Qi 

i=k L 

• ^RAxi + RBuj + Gzi) + uf Q2Ui + 

+ A^ +1 (-x i+1 + #Xj + Tuj + ©zi) 

Obviously, the last equation is similar to (3.5) if I replace RA by A, RB by 
B, G by G and all ©, Zi by ©, zj. This will lead to modified matrices Sy: 

, , k+N-l 

L = -x£ +N Sx k+N + - jP [( *7 u 7 *i)‘ 

“ i~k 

f A t R t QiRA A t R t QxRB A t R t QiG \ / \ 

B t R t QiRA B t R t QxRB + Q 2 B t R t QiG Ui + 
V G t QiRA G t QiRB G t QiG / \ / 

+ A i+i (-xi+i + $x; + Tui + ©Zi)] 

L = -x k+N Sx k+N + 

, k+N- 1 ( S n S 12 S 13 \ / X i 

+ 2 £ ( *?* n ? ) Sal S22 S23 Ui 

2 \ S 31 S 32 S 33 / \ Zi 

+ a£. x (-x i+ x + ^Xi + Tui + ©Zi)] 

Again: Qx = and Q 2 = Qj, thus Sy = Sji Vi,j S T = S 
Qi, Q 2 and S are real symmetric weight matrices. 

Keeping these substitutions in mind I can benefit from all results of the 
previous chapter. Thus, I will only state the results for the sake of com- 
pleteness: 

Pi = [ 1 - (^ T Pi + ir + Si 2 )M inv r T ^- T ] _ ' • 

• [* T P i+ i (* - rM inv (r T *- T Su - S 21 )) + Sn - 
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- S 12 M lnv (r T *- T Sn - Sai)] (4.5) 

for i = k + l...k + N — l 

Pi = [i - (* T p, +1 r + s 12 ) M lnv r T *- T ] _1 [* T Pi+1 + (* T P i+1 e + 

+ s 13 - (^ T P i+ ir + s 12 ) M inv (r T #- T Si 3 - s 23 ))zi] (4.6) 

fori = A; + l...fc + iV-“l 

u k = - (s 22 + r T P k+1 r) _1 [(s 21 + r T p k+x *) x k + 

+ (s 23 + r T P k+1 ©) z k + r T Pk+1 ] (4.7) 

In order to evaluate how good my goal is achieved depending on the sensor 
positions, I will use a set of cost functions like in the previous chapter, 
however, with the “ride comfort” as defined for this approach: 

j\ = f w2 T (t) RC-w2 T (t)dt ride comfort (4.8) 

Jt 0 

4.2 Necessary assumptions 

Even though I can “recycle” the whole calculation of the previous chap- 
ter by using the above mentioned substitutions, I have to make sure that 
the necessary matrix inversions are possible within the new system as well. 
Again, I can benefit from the preceding investigation, but have to consider 
the modified definitions: 

S 22 


S 22 = B t R t QiRB + q 2 

If (B T R T ChRB + Q 2 ) is a positive definite matrix then its determi- 
nant is unequal to zero. 

Similar to the previous argumentation for S 22 , if the terms of Qi that 
refer to nonzero parts of RB are at least positive semidefinite then 
B t R t QiRB will be, too. Again, I will assume that Q 2 is positive 
definite. Sj 2 exists 


(S^S^r-Saa) 

The transformations are similar to the previous case although there 
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are some terms of R and R T involved. The two alternative conditions 
finally read: 


B t R t Q 1 R^" 1 B + q 2 
I* + RBQ2 1 b t r t q 1 


± 0 
* 0 


[1 - (<& T P i+1 r + s 12 ) M inv r T *- T ] 

The argumentation is equal to the previous chapter since there are no 
decompositions of Sy and no © involved. 

(s 22 + r T p k+ 1 r) 

The argumentation is equal to the previous chapter since there are no 
decompositions of Sy and no © involved. 


Chapter 5 

Computational Performance 
Requirements 


At first sight, the various matrix inversions seem to be a serious problem 
when it comes to implementing the optimal control law on a DSP. This is, 
however, not the case, as will be shown in this chapter. 


5.1 Minimization of accelerations 


Equations (3.16) and (3.17) as well as the included substitution Mj nv do 
not contain the variables xj, u\ or Zj. As a result, the whole set of matrices 
P k+N • • ■ Pk+i can be computed offline and saved. 

Then, (3.19) can be simplified as follows: 


where 

Hi = 

E2i = 


p i = Fl i (^ T p i+1 + F2 i z i ) (5.1) 

i - (* T p i+1 r + Sia) M inv r T *- T ] _1 
* T p i+ i© + s 13 - ($ T P i+1 r + Sia) Mjnv (r T *- T s 13 - s 23 ) 


for i = k + 1 . . . k + N - 1 

Obviously, FI and F2 only depend on constant matrices and the known 
P k+N . . . P k +i so that those matrices can be pre-computed and stored as 
well. Furthermore, the fact that p k +N and the unknown gust z k+n . . z k+N _! 
are set to zero allows another “short cut” for the calculation: pi will (start- 
ing the recursion with i = k + N — 1) be zero as long as z; = 0, that is for 
i = k + n — k + N — 1. To save time, for 
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• 77 = 1 set directly 

Pk+i = o 

• n > 1 start the calculation with 

Pk+n-l = Flk+n-1 I^k-fn— 1 z k-fn— 1 (5-2) 

The optimal control law (3.21) can be simplified as well: 

u k = F3 (F4 x k + F5 z k + r T p k+ i) (5.3) 

F3 = - (s 22 + r T p k+1 r) _1 

F4 = S 2 i +r T p k+ i^ 

F5 = s 23 + r T p k+1 © 

Result: FI . . . F5 and P k+ N • • • P k +i can be calculated offline (before the 
experiment) and stored in the controller. During the experiment only pi, i = 
k + n — 1 ... k + 1 (or none for n = 1) need to be calculated online using 
simple and fast additions and multiplications. Slow matrix inversions are 
not required at that time. u k is finally derived using such easy operation as 
well. 


5.2 Minimization of jerkiness 


Also for this approach, offline calculations for P k +N • - • Pk+i are possible 
(equation (4.5)). The difference concerning the online part is (apart from 
the modified definitions of Sy) the fact that now © and zj come into play. 

Fii = [i- (^ T p i+1 r + s 12 )M inv r T ^- T ] _1 

F2j = * T P i+ i© + s 13 - (* T P i+ ir + Sia) M inv (r T *- T S 13 - s 23 ) 


As a result, I have to be very careful with the above mentioned “short cut” 


for p k+ i, as: 


Pi = Flj 


^ T p i+1 + E2i 



(5.4) 


The fact that p k+ N and the unknown gust z k+n . . .z k+ N-i are set to zero 
requires to start the downward calculation for i = k + n, as opposed to the 
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previous case. For n = 1, Pk+i cannot be set to zero any more. The other 
matrices FI for this approach read 

F3 = - (s 2 2 + r T P k+1 r) _1 

F4 = S 2 i + r T p k+ i* 

RS = s 23 + r T p k+ i0 
leading to the simplified control law: 

Uk = F3 (^F4 Xk + F5 Zk + r T pk+i) (5.5) 

5.3 Possible limitations 

Depending on the algorithm for matrix multiplications the computation time 
will increase with 0(q 3 ) in the worst case, where q is the dimension of 
two multiplied (square) matrices. For system models of very high order, 
the multiplication could become a problem. As shown above, the number 
of steps for the recursive calculation of Pk+i depends on n, that is the 
“distance” between the sensor(s) and the wing. For measurement far ahead 
of the airplane, the recursion could require too much time and would not be 
applicable any more. 



Chapter 6 

Implementation with 
Matlab / Simulink 


Now, I would like to implement the above-derived optimal control algorithms 
with Matlab/Simulink, which is well-established in the field of control en- 
gineering. These “experiments” will evaluate how good my goal of gust 
alleviation is achieved depending on the parameters involved. I will split 
the w r hole model in several subsystems w r hich are described hereafter. 



Figure 6.1: The test configuration and its subsystems 
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6.1 Dryden model of gust 


To receive meaningful simulation results, a realistic gust signal needs to be 
generated. I will use the Dryden model which consists of a white noise 
generator and the following second order lowpass filter: 


G(s) 


1 + \/3 UJg 
(1 + LOgS) 2 


In chapter 3 I pointed out that the gust signal should be equal within 



White Noise Dryden gust model 


Figure 6.2: The gust source using the Dryden model 

the experiments to ensure a maximum of comparability between the results. 
I proposed to sample and save a random gust function and recall it for 
every test. The noise generator of Matlab/Simulink does provide the same 
sequence for each simulation, so that a “gust data file” is not necessary. 

6.2 Airplane state space model 

Despite the fact that the controller design is based on a discrete-time model 
(and also a discrete-time cost function), the “real” system of course remains 
continuous. Both representations lead to the same states for t = kT. The 
continuous- time system, however, also provides information between these 
points in time and should therefore be used for evaluation. 

6.3 Controller 

The controller is realized for discrete-time. For this reason, all input signals, 
that is the gust right at the sensor and the current states, are sampled. Of 
course, the control algorithm itself could be realized by a m-file that is run for 
each time step t k like a program on a DSP. For the application on Simulink, 
a transfer function, however, fits better into the signal path design. As a 
side effect, the transfer function implementation allows to have an easier 
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Figure 6.3: The continuous- time model of the airplane 

insight into the controller. I will derive this for the approach of chapter 3 
first and then discuss the necessary modifications for the algorithm of 4. 

6.3.1 Minimization of accelerations 

For the case of n > 1, (5.1) and (5.2) lead to a set of coupled equations: 

Pk+n-l = Flk+n-1 F2 k+n _x z k+n-l 

Pk+n— 2 = ITk+n— 2 (^"^Pk+n— 1 d" I^k+n— 2 Z k _|_ n _ 2 ^ 

Pk+n-3 = Flk+n -3 (# T p k+n _ 2 + F2 k+n _3 Z k+n _ 3 ) 

Pk+2 = Flk+2 (^ T Pk+3 + Rk+2 z k+2) 

Pk+1 = ITk+1 (^ T Pk+2 + F2k+1 z k+l) 

Eleminating p k + n -i • • Pk+2 I obtain: 

Pk+i = Flk+ 1 * T - Fl k+2 ^ T FL k+n _ 2 $ T - Fl k+n _x- E2 k+n _i- z k+n _x + 

Fl k+ 1* T Fl k+2 4> T - ••• n k+n _ 2 F2 k+n _ 2 z k+n-2+ 


Fl k+ i^ T - Fl k+2 - 
FL k+ r 


z k+2 + 
z k+l 


F2 k+2 - 

F2k+1 - 


(6.1) 
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Now, I apply a z-transformation on this equation. To avoid confusion I will 
rename the variables first: 


Sk z k 
Pk := Pk+1 


Using the relation 

Z{ gk+i} = z l G x (z ) 


Note: G z (z) has nothing to do with the matrix G of the state space model. 
The z-transformation of (6.1) can be written as follows: 


PzW 

Pz(2) 


n— 1 

£ 

i— 1 t iiy—\ 


n n 


k+i^ 


^” T E2 k+i z l G z (z) 


n— 1 ( 

n ptk + „* T 

*- T E2 k+i z4 

i= 1 l 


J 


& x (z) 


Now, I have a filter PF z (z) with g k = z k (that is the gust directly at the 
wing) as the input and p k = p k +i as the output. The controller, however, 
will receive the signal m k of the gust measurement device, that is r? — 1 steps 
ahead: 


gk — m k _( n _x) 

^ z-transformation 

G z (z) = z _(n_1) M z (z) 


Finally the filter for m k providing Pk+i reads: 


PF z (2) 


EL'i 1 { [nUi Flk + ^ T ] ^~ T K2 k+i z’} 


Note: For the case of n — 1, p k+ i = 0 so that this filter has to be set to 
zero. The z-transformation of the optimal control law in the form of (5.3) 
with these substitutions reads: 


= F3 
= F3 


F4 X z (z) + F5 G Z (z) + r T P z (z)] 

F4 X x (z) + F5 z-( n ~^M x (z) + r T PF Z (z) M z (z)] 


U z (z) 
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= F3 F4X Z ( 2 ) + (lb2- (n - 1) + r T IT' z (z))M z (0) 

^ ' 

GF,(J) 

U x (z) = F3 [F4 X z (z) + G3F z (z) M z (z)J 


(6.2) 


G3F z (z) is a kind of “gust filter” and has the effect of a dynamic feed-forward 
part of the controller. 


GF z (z) 


r T gg { [nUi ^k + ,^ T ] *- T re k+i **} + re 


Nun n -iz n ~ 1 + Nun n -2^ n ~ 2 j hNun2^ 2 + + F5 

—— 

Nurii are matrices of the same size as G. The discrete transfer function 
module of Simulink which I would like to use is restricted to scalar input. 
Therefore, z k and have to be a scalar within this evaluation. 



Figure 6.4: The controller including the feed-forward “gust filter” for the 
case of chapter 3 


6.3.2 Minimization of jerkiness 

The following calculation for the gust part of the controller is similar to the 
case above, but I have to consider Zj. Equation (5.4) alkws to describe Pk+i 
as a linear combination of gust data: 

Pk+l = Flk+l^ T - FL k+ 2 * T ••• Flk+n— l^ T ‘ Flk+iT F2k+n Z k + n + 

Flk+1^ T - Flk + 2* T Flk+n— 1 ’ F2k+n-l’ Zk+n-l + 


n k+ i^ T - Fik+2 
Flk+1- 


Zk+2 + 
Zk+1 


E2k+2 

F2k+i- 


(6.3) 
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I have to consider that only z k . . . z k+n _i are available and z k+n . . . z k+ N_i 
are set to zero as described at the beginning. Therefore I have to be careful 
with z k+n : 

Zjc+n = 

Again, I will rename the variables and apply a z-transformation on this 
equation. 


Zk+n 


0 

Zk+n-1 


Zk+n— 1 


gk 

gk 

Pk 


= Zk 

= z k = | Z k z k _l 
= pk+l 


= Z k Z k _l J 


Using the relation 


^ {gk+i} = 
the z-transformation of (6.3) leads to: 


z l I 

z*- l I 


G z (z) 


r n - 1 


p z ( z ) — \ y. 


i=\ L \i/= 1 


n FI k+ ^ T ] ^- T E2 k+i [ z i I z*- 1 1 ] ' 


+ 


+ m U'X 


<i/= i 


B [o z"- 1 !] 71 ! 


G z (z) 


P z(z) = PF x (z)G 2 (z) 

Finally the filter for m k providing p k +i reads: 

(nt=i Flk+^ T ) ^" T F2 k+i [ Z’l z*-'l ]' 


E n— 1 

i= 1 


FF Z (z) = 


yTl — 1 


+ 




(W=i FW* T ) ^~ T E2 k+n [ O z"- 1 ! ] ' 


y n— 1 


Using these substitutions, the optimal control law in the form of (5.5) reads: 


U,(z) = F3 


= F3 


F4X z (z) + F5 


I 

z" 1 ! 


g z ( z ) + r T p z (z) 


F4 X z (z) + ^B5 z l Xl z~( n ~V + r T PF z (z)^J M z (z) 


GF.CO 
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U 2 (z) = F3[F4X 2 (z) + GF z (z)M 2 (z)] (6.4) 


Obviously the controller layout is similar to the one dealing with accelera- 
tions. 



The restriction of scalar and also applies to this case. 


6.3*3 Bounds of the controller output 

With respect to the application of my controller, it is important to know the 
upper bound of outputs. Real controllers and actuators can only operate 
within a limited range, so that a non-linearity comes into play if outputs 
exceed these limits. This would certainly reduce the performance w r ith re- 
spect to my goal. The chief problem, however, is that such non-linearities 
can jeopardize stability, w T hich is often only verified for the linear case. In 
course of the following I will determine the maximum controller output for 
a somehow “worst case”. I will assume that the initial value of the states is 
zero, as for a stable linear system this part (superposition) will decrease to 
zero anyway. Thus, the states are a linear function of gust data, so that the 
upper bound for the controller output is proportional to a kind of maximum 
of gust. To find this upper bound for all possible sequences of gust I will 
make use of norms theory, especially of H ^ which is known from robust 
controls. 

As for the definition of the used vector norms I am quite free, as long as 
some constraints are satisfied. The H 0 c norm for the discrete-time transfer 
function matrix is defined as 

F = supafF(e^)) 

oo u; > ' 

a is the greatest singular value of F*F. This norm is not trivial to calculate. 
Therefore I will make use of the command nonnhinf (sys) provided by the 
robust control toolbox of Matlab. For this reason I have to transform the 
whole controlled system, especially the “gust filter” of the controller, to state 
space. Finally I will be able to formulate the problem as 

U z (z) = F Ui (z)M 2 (z) =>■ ||u|| < F u ^ • [|m|| 
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Xz(2) = F Xi (2)M z (z) 


lix < 


m 


I will also determine the maximum for the states of the system which can 
be quite useful in view of structural displacements. 

Taking the control law in the form of (6.2) or (6.4) as my starting point: 

U Z (z) = F3 [F4X z (z) + GF z (z) M z (z)] 

= F3F4X z (z) + F3GF z (z)M z (z) 

= KX z (z) + F3GF_num z (z) z - ( n_1 )M z (z) 

= KX z (z) + F3 GF_num z (z) G z (z) 

n— 1 

U,(z) = KX z (z) + Bi^G z (z) 

i=— 1 

B; = F3 GFjmmi 

where GF_numj are the coefficients of the numerator polynomial of GF z (z) 
for the denominator being z n_1 . The inverse z-transformation leads to 


n— 1 

Uk = Kx k + Bjgk+i 

i— — l 


(6.5) 


As mentioned above, each of the “systems” F Uz (z) and F Xz (z) has to be 
described in state space, which requires to introduce new states for all known 
gust data. The corresponding dynamic matrix is a kind of shift register. 


^gk+i 


L gk4 1 


= *g x gk 


B g m k 

O I o o 
O O I o 

O O O I 


o 

o 

o 


o o o o o o 


A22 



g k -l 


o ' 


gk 


o 


gk+1 

+ 

o 


gk+n-1 _ 


I 


m k 


Now, I will transform the system equation using (6.5): 

x k+1 = ^x k + Tu k + ©g k 

71 — 1 

= 4>x k + TKx k + T Big k+ i + ©g k 

<=-i 

x k+ i = ($ + TK)x k + [ rB_! (rB 0 + ©) TBx 

? ' 

Al1 a 12 


tb,,. 


^gk+i 
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The resulting fictitious system having m k as input reads: 

+ 


Xk+X 


An A12 


x k 

Xgk + l . 

y \ 

, 

CS 

< 

o 

J 

/V 

x gk 

J 

x “»k+l 


A a it 


x «U k 



m k 


( 6 . 6 ) 


The output “y” of the system will be the vector for which I want to determine 
the maximum norm. 


for Xk 


y«u k 

:= x k = 

D all x 

= O 


= [i o] 


X a ll k 


Call* 


• for ilk 


yali k 

yall k 


D all u 


Uk 


= Kxk + B_i Bq Bi - - • B n _i j x gk 


K B_x B 0 Bi 


B n . 


X a ll k 


o 


These matrices allow to make use of a very efficient function of the robust 
control toolbox: 


sysx-ss(Aall ,Ball ,Callx>Dallx ,T) ; 
sysu-ss (Aall ,Ball ,Callu,Dallu,T) ; 
normx-normhinf (sysx) 
nonnu*normhinf (sysu) 


The whole m-file ( Hinf.calc.m ) is included in the appendix. Finally I obtain 
for the boundaries of x and u: 


Ml < 

ii*ii < 



IMI 

l|m|| 


where 



— normu and 

oc 



= normx. 

oc 
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Now, the question is how good these “upper bounds” of the above stated 
inequalities are with respect to my application. Therefore I would like to 
review the results for a general case: 


Y z (z) = F z (,z)U z (z) => ||y |I < 

I have to assume that my input signal is band limited: 

TJ (u^) — 0 V u; > 'jJfxuxx 


u 


I would like to give the definition of the H-infinity norm again: 


= sup o (*(e , ’“')) = a (F(e ja, '"* p ) 


In case jj sup is greater than u} max , the inequality for ||y|| is a very rough 
estimate and a constant c exists such that 


llyll < 

C < 


C ■ Hull 

11*11 


Thus || F ||oc is not the supremum of c for this particular input signal. 
Obviously, it is necessary to regard the results provided by Hinf-calc.m in 
view of the application. If the eigenfrequencies of the controlled system are 
within the bandwidth of the input signal (gust) the H-infinity norms are a 
good estimation. 


6*3*4 Performance indices by means of H-infinity 

With respect to evaluation of the controller designs it is also interesting to 
compare similar H D 0 norms for acceleration and jerk, respectively, of the 
controlled system and the system without the controller. These norms can 
also help to guarantee that specifications e.g. for maximum acceleration are 
observed. 


An,(z) 

= F ACi (z)M z (z) 


Kll 

< 

||Fan||oo 

• ||rn|| 

no contr. 

A Ci (2) 

= FAru(z) M z (z) 


IM 

< 

||F a c||oo ' 

• Ml 

controlled 

Jn,(z) 

= F JCi (z)M z (z) 


II Jn 1 ! 

< 

IIFjnlloc ' 

IHI 

no contr. 

Jc(*) 

= F Jni (z)M z (z) 


IIJcll 

< 

IIFjclloc • 

Mil 

controlled 


These calculations are also part of Hinf-calc.m and use the fictitious system 
described by A a n, B a u and x a n k for the controlled system norms. As for 
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the system without a controller, I replace A a n by A a u. I will set a = Lx for 
the first and j = L w2 for the second approach. L selects the states related 
to the fuselage as this is in line with my design goal, ||L|] = 1 is assumed. 


x gk+1 = x gk + Bgirik as for the controlled case 


A22 

Xk+i = $x k + ®Sk 

x k+l = ^ x k + 0 © O ] 

- v “ V 

An 


Xgk 


+ i 


A12 

Therefore, for the system without a controller I obtain: 


x k +i 


>< 

I-* 

M 

1 >« 
H* 

K5 

X k 

x gk+l 

V V 

S V 

0 a 22 J 

S- Xgk V 



m k 


(6.7) 


Xall k+1 A aJ i X * ,l k ^all 

Now, I can determine the various C and D matrices so that the output of the 
fictitious transfer function equals the variable for which I want to calculate 
an upper bound. 

• accelerations of the system without a controller 


a n k = Lx k 

= L (Ax k + Gg k ) 
yaii k = [la 0 LG O ] x a ii k 

' V ' 

C a ll an 

L) a n an O 

Use Aaj| for the dynamic matrix. 


• accelerations of the system with the controller 


ack 


Yallk 

D all ac 


Lx k 

L (Ax k + Bu k + Gg k ) 

LAx k + LB ^Kx k + [ B_i Bq Bj • B n _i j x gk ^ + LGg k 


(A + BK) BB_i (BB 0 + G) BBi 


c alla 


BB n . 


x all k 


o 


Use Aan for the dynamic matrix. 
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jerkiness of the system without a controller 

j„ k = L w2 k 

. r ~ 1 x k 


W2 = -M2 1 [ 
w2 = -M2" 1 [ 


K2 D2 O ] 

K2 D2 O ] x k 


x k 


— L M2 1 

K2 D2 O 

X k 

A 0 G O ] x a n k 


— LM2 -1 

K2 D2 O 

f A 0 G O ] 


L J 

l J 


C<lll Jn 

D a n Jn = O 

Use A-aii for the dynamic matrix. 

• jerkiness of the system with the controller 

j Ck = L w2 k 
w2 k = Rx k + G2g k 

~ RAxj( + RBu k + G 


gk 

gk-X 


R ( A -f BK) x k + RB B_i Bo B^ 


B n -i x g k + 


+ G 


I 

_ 

O 

O 

gk + G 

I 


gk— 1 


\v2 k 


R(A + BK) (rBB_! + g[o I ] T j (rBB 0 + g[i O 
RBB n _i 


RBB 


x all k 


R(A + BK) (RBB-i + G 


yallu = L 


RBBi • • • RBB n _j x a n k 

y a ll k = C a ll Jc X a l| k 

D a ll Jc = O 

Use A a ii for the dynamic matrix. 


r i t \ / 

[ ° 1 ] ) ( 


RBB 0 + G 


I O 
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The difference between ||Fan||oo and ||F ac ||cx> as well as ||Fj n ||oc and HFjclIsc, 
respectively, provides information about the beneficial effect of the particular 
controller for the worst case of gust. I will divide these differences by the 
norm for the controlled case to obtain an index telling how many percent 
the system behavior is worse without the controller: 


improvement (a) == 


improvement (j) = 


F an 

oc 

F ac 

OO 


F ac 


OO 


• 100 % 



• 100 % 


In the previous section, I discussed the meaningfulness of the boundaries 
calculated using H-infinity norms. Since the performance indices are based 
on the same theory, I should consider the bandwidth of the input signal for 
the particular application as well. 


6.4 Evaluation module 

This subsystem contains the cost functions J 3 , J 2 and J\ or J\ for each test 
to provide easy to compare results. The scopes show the increase of cost 
during the experiment. Hereafter is shown the subsystem for the approach 
of chapter 3 and chapter 4. To calculate the jerkiness and the ride comfort 
within the second approach I make use of a derivative element . I would like 
to point out that such an element is not used for the controller itself. 
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Figure 6.5: The evaluation subsystem for acceleration control 



Jerkiness of w2 



Figure 6.6: The evaluation subsystem for jerkiness control 
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6.5 M-files 

All constants within the evaluation are defined in const.def.rn, that is: The 
sytem (described by A, B, G), the weight matrices (S, Qi, Q2), the sam- 
pling period T, the number of available gust data n, the length of the finite 
horizon N and the parameter uj 9 for the Dryden model of gust. 
const-calc-l.m, consLcalc.2.m respectively, calculate all variables that de- 
pend on the definitions in const-def.m. The calculation of the various H 0 0 
norms is implemented in Hinf.calc.m. Finally, pzmap.calc.m plots the pole- 
zero map for the system with and without the controller to verify stability. 
All of these files are included in the appendix. 



Chapter 7 

Preliminary Simulation 


7.1 Test configuration, system and actuator 

To examine the performance of the two different controllers and their im- 
plementation a very simple oscillator (pendulum) will be subject to control. 
Its differential equation reads: 


w + -w + W = J 
z 

The step response of the system without a controller is shown in figure 7.1. 
The disturbance (gust) and the actuator lead to a force / that has influence 
on the oscillation. Often there is a delay between the occurrence of a distur- 
bance and its effect (here the unwanted oscillation), as well as between the 
“command” of the controller and the moment the countermeasures really 
come into effect (because of actuator delays). With regard to this problem 
information of future disturbance should be able to increase the controller 
performance. To stress this effect, I chose a slow PT2 actuator. Its transfer 
function reads: 

G ^ = 0.04s 2 + 0.2s + 1 

The step response is shown in figure 7.2. I also tried a PT1 actuator for the 
controller dealing with jerkiness, since the difference of order is smaller in 
that case: 

G < s > - 7TT 

The step response is shown in figure 7.3. 
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Step Response 

It.'; 



^0 ’S at' a-: 

Time (sec.) 


Figure 7.1: Step response of the pendulum itself 


Step Response 



Time (sec.) 


Figure 7.2: Step response of the PT2 actuator 
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Step Rwpon** 


f JJ> V; 



Tim® (»®c.) 


Figure 7.3: Step response of the PT1 actuator 
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7.2 Minimization of accelerations 

Since this approach is not critical regarding the disturbance g k (or z k ) I 
can use a simple step function that allows easily comprehensible results to 
determine whether the control algorithm works as expected. 

The main criteria for the comparison of the cases n = 1 and n = 30 is of 
course the ride comfort index at the end of the simulation. I tried to keep 
the magnitude of the controller outputs within the same range, so that the 
performance does not depend on that. I also calculated the control energy 
for the sake of completeness, although the energy will be - as opposed to the 
output range - not be a problem to most applications. The displacements 
should of course not become too large, even though it is not the purpose of 
this controller to keep them as small as possible. Obviously, the results are 
much better, if more gust data are available. Then, as can be seen from fig- 
ure 7.6, right graph, the controller begins its countermeasures already before 
the disturbance encounters the system, what leads to lower accelerations. 


Hid* Contort Fid* Contort 




Figure 7.4: “Ride comfort” as the integral of squared accelerations for n=l 
and n=30 








Figure 7.6: Controller output for n=l and n=30 
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AcctfmAon «T Acr^nTton 



Figure 7.7: Acceleration of the pendulum for n=l and n=30 


MqUoamam w D«pt*c*m*** 



Figure 7.8: Displacement of the pendulum for n=l and n=30 
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Figure 7.9: “Aggregated displacements 1 ’ as the integral of squared displace- 
ments of the pendulum for n=l and n=30 




Figure 7.10: Control energy for n=l and n=30 
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7.3 Minimization of jerkiness 

Dealing with the jerkiness is critical with respect to the time-derivative of the 
disturbance. It is important that the disturbance function is differentiable 
for every time. Hence, I choose the first half wave of a sin 2 (u;t) function, 
which also has the advantage, that its derivative is zero for t = 0. 

7.3.1 PT1 actuator 


FM* Comfcrl 



Figure 7.11: “Ride comfort” as the integral of squared jerkiness for n=l and 
n=30 
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<M m, t Qutt * 




Figure 7.12: Measurement data and disturbance for n=l and n=30 



Figure 7.13: Controller output for n=l and n=30 







Figure 7.15: Acceleration of the pendulum for n=l and n=30 









Figure 7.16: Displacement of the pendulum for n=l and n=30 



Figure 7.17: “ Aggregated displacements” as the integral of squared displace- 
ments of the pendulum for n=l and n=30 
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Figure 7.18: Control energy for n=l and n=30 
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7.3.2 PT2 actuator 

For the jerkiness control with a second order actuator I have to assume at 
least n = 2. The gust filter calculated for this approach (PT1 as well as 
PT2) is a weighted sum of gust data m k _ l7 m k , . . . , m k+n _i where t k is the 
current point in time. As opposed to the case of a first order actuator, here 
the factorial for m k _i is always equal to zero, since B T R T = O. In view of 
the approximation for g k I need at least two consecutive gust data. Thus 
n = 1 is not applicable. Also for this simulation, an improvement is visible 
when increasing the available data. 


Ada Cemtart RhACW*r1 



Figure 7.19: “Ride comfort” as the integral of squared jerkiness for n=2 and 
n=30 
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Figure 7.21: Controller output for n=2 and n=30 










Chapter 7. Preliminary Simulation 


6 



Figure 7.25: “Aggregated displacements’" as the integral of squared displace- 
ments of the pendulum for n=2 and n=30 
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Figure 7.26: Control energy for n=2 and n=30 






Chapter 8 

Application on a 
Wing-Fuselage System 

8.1 Test configuration, system and actuator 

For the main experiment I choose a two degree of freedom system. Of course, 
with respect to my theory, I could use a more difficult model since neither 
the number of actuators, nor the degree of freedom is limited. However, 
this 11 simplified airplane” consisting of a wing, a fuselage and one first order 
actuator is sufficient for evaluation, as I am interested in qualitative rather 
than in exact results. The sketch of figure 8.1 illustrates the variables and 
constants of the state space representation: 


m 0 w 2 j D\ -D\ w< 2 ] 

0 M u>2 2 —D\ D\ + Di u>2 2 



K 2 w 2 B 2 wi 


For the determination of the spring and damping constants, I made the 
following assumptions: 

• The mass of the fuselage is 100 times higher than that of the wing 
which I set to 1. 

• The natural frequency (mode) of the wing itself is about 4Hz and the 
damping is 10%. 
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Figure 8.1: Sketch of the wing-fuselage system 


• The fuselage itself has a mode at about 1Hz and a damping of 2.5%. 

This leads to the following set of constants: 

Mass m = 1 M = 100 
Spring K\ = 700 K 2 = 4000 
Damping D\ =20 D 2 = 0.3 

The natural frequency of the fuselage is clearly visible at the bode plot (figure 
8.2). For the step response of the coupled system (figure 8.3) the system 
input is an external force (the gust) to the wing. 

The actuator model contains all dynamic parts including the servo-motor 
and the mechanism moving the ailerons or flaps. Therefore the time constant 
is higher than for a typical servo itself. The step response is shown in figure 
8.4. 

I will refer to the Dryden model of gust, as described above. Figure 8.5 
sIkws the Bode plot for the case of u; gust = 0.1. 




Amplilude Pha»e (deg); MagnitucJe (dB) 
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Bode Diagrams 



Frequency (ratWsec) 


Figure 8.2: Bode plot of the wing-fuselage system 


Stop Response 

v-icr tin, 
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Figure 8.3: Step response for the wing-fuselage system 
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Step Response 



Time (aec } 


Figure 8.4: Step response of the first order actuator 


Bode Diagram* 

'■'icr 



Frequency (rad's ec) 


Figure 8.5: Bode plot of the filter for the Dryden model of gust 
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In order to know the system behavior without the controller, I would 
like to give the same graphs that will be discussed in the following two 
experiments also for the system without a controller. Of course, this is 
not very realistic, as there will always be a controller - even though its 
purpose is just to reduce the displacements - so that the damping would be 
higher. However, I think it contributes to the understanding of the system. 
In anticipation of the main experiments, I would like to give the norms for 
acceleration and jerkiness for this system without a controller: 

iiFalloo « 179 

IlFjlloo ~ 1129 

The eigenfrequencies of the system are around 1Hz, 4Hz respectively, and 
thus within the band width of the Dryden model of gust (see 8.5). Although 
||g|| is approximately 8.5 for this simulation, the acceleration and the jerk- 
iness are much smaller than the limit given by means of the corresponding 
H-infinity norm. On the one hand, the applied “turbulence” was certainly 
not the worst case, which would be a sinus wave matching a system mode. 
On the other, the simulation, especially the period of disturbance, was per- 
haps not long enough: The displacements (see figure 8.8, left graph) seem 
to be on the increase even at the end of the disturbance, what suggests that 
a kind of “steady state” on a higher level is not yet reached. Thus, for 
quantitative comparisons with the following simulations I rather refer to the 
H-infinity norms. For the evaluation among different sensor positions the 
test period will prove suitable. 
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Figure 8.8: Displacement of the fuselage and the gust applied to the wing 
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8.2 Minimization of accelerations 

First of all I would like to refer to the pole- zero maps (figures 8.9 and 8.10) 
to check for stability. For both values of n the eigenvalues are within the unit 
circle, the controlled systems are stable. I don’t have to worry about the 
feed forward part of the controller, since for a stable linear system, bounded 
“disturbances” - the gust and the feed forward component which is only a 
weighted sum of a finite number of bounded gust data - can never lead to 
unbounded states and jeopardize (bibo) stability. 


Pole-zero map 



Real Axle 


Figure 8.9: Pole-zero map: system without (green) and with the controller 
for n=l (blue) 

Now, I would like to discuss the H- infinity norms as calculated by Hinf.calc.m: 


Index 

n=l 

n=30 

ilfaclloc 

0.5070 

0.1737 

1 1 F an 1 1 oc 

179.6118 

179.2839 

IIFulloc 

32.2925 

32.3339 

IlFxlloc 

10.2074 

10.0359 

Improvement a 

353.2302 

1031.2 
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Pole-zero nwp 



Figure 8.10: Pole-zero map: system without (green) and with the controller 
for n=30 (blue) 


The values of ||Fan||oc for n = 1 and n = 30 should be equal, the 
difference is due to inaccuracies of the calculation. For both, the system 
with and without a controller, the modes are around 1Hz, 4Hz respectively, 
and thus within the band width of gust. I would like to give the following 
upper bounds for ||g|[ « 8.5 with respect to the inequalities stated in section 
6.3.3: 

bound for n=l n=30 

IJadl 4.31 1.48 

||u|[ 274.49 274.84 

The controller output remains clearly below these limits within my sim- 
ulation, but for a more critical gust function some peaks could be higher. 
The accelerations for n — 30 are already close to the theoretical borderline 
for this particular controller but for the case of n = 1 the performance will 
deteriorate in the worst case making the advantage of additional gust data 
even more obvious. 

The main criterion for my evaluation is the “ride comfort” which I de- 
fined (for this approach) as the integral of the squared accelerations. Ev- 
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ery simulation lasts 25 seconds and I apply the same gust sequence for 
1 $ < t < 11 s to the system, so that the results are easily comparable. 

It is unquestionable that already the controller for n = 1 leads to a signifi- 
cant improvement of all considered criteria. Thus, I will concentrate on the 
question whether the additional information available in the case of n = 30, 
which equals data almost about the next 1.5 seconds (T=0.05s), leads to a 
remarkable increase in performance. 

I will give the results of this comparison in tabular form. I w T ould like to 
remind of the fact that a comfortable ride has a low (!) ride comfort index. 
Increase/Decrease of the results for n=30 compared with the case of n=l: 


Ride comfort index -58% 

Controller Output -1% 

Accelerations -24% 

Displacements -47% 


Aggregated displacements -69% 

Note: The percentages of the table for the accelerations, the displacements 
and the controller output refer to the maximum absolute value of each graph. 


FUm Canton ffcto Canton 



Figure 8.11: “Ride comfort” as the integral of squared accelerations for n=l 
and n=30 
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Figure 8.16: “Aggregated displacements” as the integral of squared displace- 
ments of the fuselage for n=l and n=30 




Figure 8.17: Control energy for n=l and n=30 
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8.3 Minimization of jerkiness 

Also for this approach, I will check for stability in the first place. As can be 
seen from the pole-zero maps (figures 8.18 and 8.19) for both values of n the 
eigenvalues are within the unit circle and the controlled systems are stable. 


Pole-zero map 



Figure 8.18: Pole-zero map: system without (green) and with the controller 
for n=l (blue) 

The H-infinity norms for my two approaches show the same trend. The 
results provided by Hinf.calc.m read: 


Index 

n=l 

o 

CO 

II 

a 

!|Fjc||oo 

4.0383 

3.2450 

llfjnlloo 

1129.9 

1128.9 

l|F„||oo 

21.3451 

19.1240 

IlFxlloc 

14.4614 

10.0953 

Improvement j 

278.7991 

346.9065 


The values of ||Fj n ||ot are slightly different due to inaccuracies of the 
calculation. The modes are off course still around 1Hz, 4Hz respectively. 
The upper bounds for j[j c |[ and ||u|| are given for [|g|| « 8.5: 
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Pole-zero map 



Figure 8.19: Pole-zero map: system without (green) and with the controller 
for n=30 (blue) 


bound for n=l n=30 1 

||j c || 34.33 27.58 

jju|| 181.43 162.55 

Again, the controller outputs remain far below these limits, but the maxi- 
mum of jerkiness within my simulation reaches about 50% for both values 
of n. 

Comparing the graphs for n — 1 with the results of the system on its 
own, there is no doubt that the controller is able to reduce jerkiness. Again, 
the question is whether the extension towards n = 30 allows a significant 
increase in performance. The results of the comparison between n=30 and 
n=l read as follows: 


Ride comfort index 

-51% 

Controller Output 

-2% 

Jerkiness 

-32% 

Accelerations 

-40% 

Displacements 

-48% 

Aggregated displacements 

-73% 
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Note: The percentages for jerkiness, accelerations, displacements and con- 
troller output within the given table refer to the maximum absolute value 
of each graph. 

I would like to recall the “Fatigue Life Correlation” diagram that I dis- 
cussed at my introduction. Although the improvement due to the additional 
gust data is not as evident as in the case of the first approach or within my 
preliminary simulations, for the jerkiness, a reduction of 32% can still lead 
to a substantial increase in life time. For the example of vibration fatigue 
of SMT solder joints it is more than a factor 100. 




Figure 8.20: “Ride comfort” as the integral of squared jerkiness for n=l and 
n=30 
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Figure 8.25: Displacement of the fuselage for n=l and n=30 



Figure 8.26: “Aggregated displacements’" as the integral of squared displace- 
ments of the fuselage for n=l and n=30 
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Figure 8.27: Control energy for n=l and n=30 
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8.4 Discussion of appropriate sensor positions 

I already gave an affirmative answer to the question whether moving the 
measurement device ahead of the wings has a remarkable beneficial effect. 
Indeed, that should be the case for almost every system and I showed the 
improvement for the particular case of a wing-fuselage system. However, 
the question remains how far ahead it is appropriate to place the sensors. 
There is no general answer to this question, and the possible improvement 
is highly dependent on the particular system as well as the actuator that is 
used. In principle, additional information should be of more benefit if the 
delays between the encounter of gust and the undesirable effects on the one 
hand and the delays between the start of countermeasures and their effect 
on the other are large. 

One possibility to decide on the sensor position is to look at the feed for- 
ward part of the controller in the form of a transfer function, as described in 
chapter 6. This representation allows an insight in the control algorithm, as 
the weights of all available gust data are shown. Starting from the calcula- 
tion for the sensor far ahead (n = 30 in this case), I can examine how many 
of the highest powers of z I could neglect since their coefficients are too small 
to have a major effect on the controller performance. The highest remaining 
power of z defines the upper bound of “reasonable 11 improvements. 

To illustrate this strategy I would like to give the feed forward transfer 
function for the system of my preliminary simulation. The approach of 
minimizing the jerkiness using a PT2 actuator (section 7.3.2) provides the 
best example for that purpose: 


Vff x (z) = F3GF z (z)M z (z) = 

O.OOI62 30 - 0.0012z 29 - 0.0041z 28 - 0.0068z 27 - ... - 0.0056z 13 - 

- 0.0007z 12 + 0.0050Z 11 + 0.01 12z 10 + 0.0167z 9 + 0.0167z 8 - 

- 0.0061z 7 - 0.1148z 6 - 0.5375z 5 - 2.101z 4 - 7.795z 3 - 


- 28.45Z 2 + 29.842 


2" 30 M z (2) 


It seems that the gust data of gk . . . gk+s have the most influence on the 
controller output, so one could try whether the performance for n = 6 is 
already sufficient. 

For the wing-fuselage system, however, the coefficients of even the high- 
est powers of 2 are significantly different from zero, what shows that the 
gust data in the distant future still help to improve ride comfort in the sense 
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of the particular approach. 


The feed forward transfer function for the experiment of minimizing ac- 
celerations reads: 


U ff .(z) = 

0.4047z 29 + 0.2655z 28 + . . . + 0.7933z 15 + 1.176z 14 + 1.476z 13 + 

+ 1.652z 12 + 1.676Z 11 + 1.534z 10 + 1.215z 9 + 0.7121z 8 + 0.0505z 7 - 


- 0.6693z 6 - 1.308z 5 - 1.854z 4 - 2.508z 3 - 3.341z 2 - 


- 3.696z - 2.491 


z" 29 M z (z) 


The corresponding function for the approach of minimizing jerkiness: 


U ff ,(z) = 

0.2141z 30 + 0.1484z 29 + 0.0572z 28 - 0.0506z 27 ± . . . - 0.2906z 6 + 

+ 0.6206z 5 + 0.5423z 4 - 3.459z 3 - 8.832z 2 - 0.8890z + 0.4897 ■ 

• z- 30 M z (z) 


Chapter 9 

Final Remarks 


9.1 Conclusions 

The simulations proved that both controller designs are able to accomplish 
their individual design goal. It turned out that the main constraint for the 
application of gust measurement ahead of the wings is the availability of 
suitable measurement devices, their accuracy and the additional expenses 
for this technology rather than a limit of improvements close to the wings. 

Future work in the field of direct gust measurement should deal with 
relaxed constraints for the accuracy of data to find an optimum between re- 
liability of data and the prediction interval. The major task, however, is the 
further development of sensor technology. Aside from improvements regard- 
ing range and accuracy, the decoupling from the vibrations of the fuselage 
or the wings - wherever the device is mounted - is essential. Otherwise the 
measured data are also a function of system states and the feed forward 
property cannot be assumed anymore. 

Even though this investigation equates the disturbance to the system 
with atmospheric turbulence, both control algorithms are applicable to any 
kind of ride comfort problem provided that the given set of constraints is 
satisfied. 
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Appendix A 

Matlab M-Files 


A.l const-def.m 

This is just an example. The m-file is adapted to each experiment. 

7. Blu-fl -> [(Mlvl” + )Dlwl’+Klvl»n] 

'/, -> wl eigther first or second order diff. eqn. e.g. actuator 

7. B2wi SE f 2 -> [ M2w2’>+ D2w2 * +K2w2«f 2] 

7, -> w2 second order diff. eqn. 

7. x* [v2 w2 5 wl wl ’ ] ~T 

7, M,D,K all Bquare matrices 
7, "actuator" 

7. M1=[J; 7. if applicable 

Dl* [0 . 2] ; 

Kl* [1] ; 

B1-E0.5] ; 7. f l=Bl*u 

7, "system" [wing fuselage] *T 
M2* [1 0; 

0 100 ] ; 

D2* [20 -20; 

-20 20.3] ; 

K2* [700 -700; 

-700 4700] ; 

B2= [1 ; 0]; 7. f2=B2*vl 

G-[0;0; 10;0; 0]; 7. gust effect on x- [w2 ; w2’ ;wl ;wl ’] ~T, 

7. x*[w2 ; w2 1 ;wl] *T repectively 

S* [1000000000 0 0 0 0; 

0 100000 000; 

00000; 

0 0 0 0 0; 7. weight matrix for final value of x: 

0 0 0 0 0]; 7. x_(k+N) must be symmetric 

Q1*[0 0000; 

00000; 

0 0 1 0 0; 7, weight matrix for x r must be symmetric 
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0 0 0 10000000 0; '/♦ and soma terms pos.def. 

0 0 0 0 0 ]; 

Q2-1000; % weight matrix for u must me symmetric 

X and pos . def . 


RC-[0 0 0 0 0; 
00000; 
00000; 
000 1000 0 ; 

0 0 0 0 0 ]; 

AD- [0 0 0 0 0; 

0 1000 000 ; 

00000 ; 

00000 ; 

0 0 0 0 0 ]; 


% only lower right term unequal to zero 


% only upper left term unequal to zero 

*/. Evaluation of the control energy size as Q2 


T-0.05; 
n-30; 
N-40 ; 


'/♦ sample period 

'/, number of available gust data 
X length of the finite horizon, N>-n 


v_gust-0 . 1 ; 


*1 parameter for Dryden gust model 


A. 2 const-calc-l.m 

Dlsize*size(Dl) ; 

D2size-size (D2) ; 

Blsize-size (Bl) ; 

if exist ( J M1 * , *var J )--l 

A= [zeros (D2size) eye(D2size (1) ) zeros(D2size(l) ,Dlsize(2)) zeros(D2size (1) ,Dlsize (2)) ; 
-inv(M2)*K2 -inv(M2)*D2 inv(M2)*B2 zeros (D2size (1) ,Dlsize(2) ) ; 

zeros (Dlsize (1) , D2size (2) ) zerosCDlsizeCl) ,D2size(2)) zeros (Dlsize) eye (Disize (1) ) ; 
zeros(Dleize (1) ,D2size(2)) zeros (Dlsize(l) ,D2size(2)) -inv(Ml)*Kl -inv(Ml)*Dl] ; 
B-[zeros(2*D2size(l)+Dlsize(l) ,Blsize(2) ) ; inv(Ml)*Bl]; 

V, case of wl: 2nd order diff. eqn. 
else 

A- [zeros (D2size) eye(D2size(l) ) zeros (D2size(l) ,Dlsize(2)> ; 

-inv (M2) *K2 -inv(M2)*D2 inv(M2)*B2; 

zeros (Dlsize (1) ,D2size(2>) zeros(Dlsize (1) ,D2size(2)) -inv(Dl)*Kl] ; 

B= [zeros (2*D2size (1) ,Blsize(2) ) ; lnv(Dl)*Bl] ; 

*/, case of wl: 1st order diff. eqn. 

end 

’/, G already defined in const_def 
Asize-size (A) ; 

Bsize=size (B) ; 

Gsize-size (G) ; 

Phi-expm(A*T) ; 

Gamma=inv (A) * (Phi-eye (Asize ( 1) ) ) *B ; 

Theta-inv (A)* (Phi-eye (Asize (1)) ) *G ; 

Sll-A 1 *Q1*A ; 

S12=A 1 *Q1*B ; 

S13=A ’ *Qi*G ; 
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S21-B'*Q1*A; 

S22«B , *Q1*B+Q2; 

S23-B ’ *Q1*G; 

531- G 1 *Q1*A ; 

532- G’*Q1*B; 

533- G ’ *Q1*G ; 

Minv=inv (Gamma* * in v (Phi ’ )*S12-S22) ; 

% now calculating P_i for i-k+N...k+l, here k sat to 0 
% P_i is of the same size as A 
P-zeros (Asize (1) , Asize(2), N) ; 

P( : , : ,N)-S; 
for l-N-i:-l:l 

PC: , : .D-invCeyeCAsizeCl) )-(Phi >*P( : , : 1 l+l)*Gamma+S12) ^MinvfrGamma’ *inv (Phi * ) ) 

* (Phi ’ *P( : , : ,1+1) *(Phi-Gamma*Minv* ( Gamma' *inv (Phi 1 ) *S11-S21) )+Sll-S12 
♦Minv* (Gamma ’ *inv(Phi * ) *S11-S21) ) ; 

end 

l now calculating Fl_i and F2_i for i=k+n-l . . .k+1 , 

% sufficient as only p_ (k+n-1) . . .p_ (k+1) calculated, 

X here k set to 0 

’/, Fl_i is of the same size as A, F2_i is of the same size as G 
if n>l 

Fl=zeros(Asize (1) , Asize(2), n-1) ; 

F2-zeros(Gsize(l) , Gsize(2), n-1); 
for l=n-l:-l:l 

Fl ( : , : ,l)-inv(eye(Asize(l))-(Phi ’*P( : , : , l+l)*Gamma+S12)*Minv*Gamma , *inv(Phi * ) ) ; 
F2(: , : , 1)-Phi * *P( : , : , l+l)*Theta+S13-(Phi * *P( : , : , 1+1) *Gamma+S12) 

♦Minv* (Gamma ’ *inv (Phi * ) *S13-S23) ; 

end 

end 

X matrices required for controller 
F3— inv (S22+Garama ’ *P ( : , : , 1) *Gamma) ; 

F4-S21+Gamma > *P( : , : ,l)*Phi; 

F5 iC S23+Gamma , *P(: , : , l)*Theta; 

GF_num=zeros (Bsize (2) ,n) ; ’/, trows like u_k, n columns for z“ (n-1) . . . z"0 

GF_den-zeros (1 ,n) ; 

GF_den(l , 1)*1 ; % denominator is z~(n-l) 

for 1-1: n-1 ’/, coefficient for z*l, for n-1 the loop is skipped ! 

FProd-eye (Asize (1) ) ; 
for r=l:l 

FProd=FProd*Fl (:,:,r) *Phi * ; 

end 

GF_num( : ,n-l)-Gamma’ *FProd*inv(Phi * )*F2( : , : ,1) ; 

X coefficient for z~l is located at column (n-1) 

end 

GF_num( : ,n)-FB ; 
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A. 3 const-calc-2.m 

Dlsize^size (Dl) ; 

D2size-size (D2) ; 

Blsize*size(Bl) ; 

if exist (‘Ml* , ‘var’)”! 

A* [zeros (D2size) eye (D2size(l) ) zeros (D2size(l) ,Dlsize(2)) zeros (D2size (1) ,Dlsize(2)) ; 
-inv(M2)*K2 -inv(M2)*D2 inv(M2)*B2 zeros (D2size(l) ,Disize(2)) ; 

zeros (Dlsize (1) ,D2size(2)) zeros(Dlsize(l) ,D2size(2) ) zeros (Dlsize) eye (Disize (1) ) ; 
zeros (Dlsize (1) ,D2size(2) ) zeros (Dlsize(l) ,D2size(2)) -inv(Ml)*Kl -inv(Ml)*Dl] ; 

B* [zeros (2*D2size(l)+Dlsize (1) ,Blsize (2)) ; inv(Ml)*Bl]; 

*/, case of wl : 2nd order diff. eqn. 
else 

A* [zeros (D2size) eye(D2size(l)) zeros (D2size (1) , Dlsize (2)) ; 

-inv(M2)*K2 -inv(M2)*D2 inv(M2)*B2; 

zeros (Dlsize (1) ,D2size(2)) zeros (Dlsize(l) ,D2size(2)) -inv(Dl)*Kl] ; 
B*[zeros(2*D2size(l) ,Blsize(2)) ; inv(Dl)*Bl]; 

*/, case of wl: 1st order diff. eqn. 

end 

’/» G already defined in const _def 
Asize*size (A) ; 

Bsize*size(B) ; 

Gsize*size (G) ; 

Phi-expra(A*T) ; 

Gamma=inv (A) * (Phi-eye (Asize (1) ) )*B ; 

Theta=inv(A)* (Phi-eye (Asize (1) ) )*G; 

R-A(D2size(l)+l :2*D2size(l) , : ) ; '/. the rows of A refering to W2’* 

Thet at* [Theta zeros(Gsize)] ; 

Gt*R* [G zeros (Gsize)] +G (D2size (1) +1 :2*D2size (1) , : )* [eye(Gsize (2)) -eye (Gsize (2) )]/T; 
Thetatsize«size(Thetat) ; 

S11*A ’ *R , « i Ql*R< t A ; 

S12*A’*R , *Q1*R*B; 

S13*A y *R * *Ql*Gt ; 

S21»B>*R , *Q1*R*A; 

S22»B , *R , *Q1*R*B+Q2; 

S23*B } *R’ *Ql*Gt ; 

S31*Gt , *Q1*R*A; 

532- Gt’*Ql*R*B; 

533- Gt’*Ql*Gt; 

Minv*inv (Gamma* *inv(Phi 1 )*S12-S22) ; 

'/, now calculating P,i for i*k+N...k+l, here k set to 0 
'/, P_i is of the same size aB A 
P-zeros (Asize (1) f Asize(2), N) ; 

P( : , : , N)=S; 
for 1*N-1:-1 : 1 

P(: , : , l)*inv (eye (Asize (1))- (Phi '*P(: , : ,l+l)*Gamma+S12)*Minv*GaI^raa , *inv(Phi , )) 

* (Phi ’ *P ( : , : ,1+1)* (Phi -Gamma*Minv* (Gamma ’*inv (Phi 1 )*S11-S21) )+Sll-S12 
*Minv* (Gamma ’ *inv (Phi * ) *S11-S2l) ) ; 

end 
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X now calculating Fl_i and F2_i for i*k+n. . .k+i , 

X sufficient as only p_(k+n) . . .p_(k+l) considered, 

'/, here k set to 0 

X Fl_i is of the same size as A, F2_i is of the same size as Thetat 
Fl*zeros (Asize (1) , Asize(2), n) ; 

F2-zeros (Thetatsize (1) ,Thetatsize (2) f n) ; 
for l-n:-i:l 

Fi(: , : ,l)*inv(eye(Asize(l))-(Phi , *P( : , : ,l+l)*Gamma+S12)*Minv*Gamma , *inv(Phi > )) ; 
F2 ( : , : ,l)-Phi**P(: , : ,l+l)*Thetat+S13-(Phi ’ *P( : , : , l+l)*Gamma+S12) 

•Minv* (Gamma 1 *inv (Phi’ )*SI3-S23) ; 

end 

X matrices required for controller 
F3«- inv(S22+Gamma , *P(: , : , l)*Gamma) ; 

F4-S2 i+Gamma ’ *P ( : , : , 1) *Phi ; 

F6»S23+Gannna * *P ( : , : , 1) *Thet at ; 

GF_num-zeros(BBize(2) ,n+l) ; X #rows like u_k, n+1 columns for z"n...z"0 
GF_den=zeros ( 1 , n+ 1 ) ; 

GF_den(l , 1)*1 ; X denominator is z~n 

FProd-eye (Asize (1) ) ; 

for l“l:n-l X coefficient for z"l 

FProd«FProd*Fl ( : , : , 1) *Phi 5 ; 

GF_num( : ,n-l) S£ GF_num( : .n-D+Gamma’ *FProd*inv(Phi * )*F2( : , : ,1) 

* [eye(Gsize(2)> ; zeros (Gsize (2) ,Gsize(2))] ; 

X coefficient for z“l is located at column (n+1-1) , 

X additional shift left because of {}z ! 

GF_num(: ,n- (1-1) )=GF_num( : ,n- (1-1) ) ♦Gamma’ *FProd*inv(Phi * ) *F2( : , : ,1) 

* [zeros (Gsize (2) , Gsize (2)) ;eye (Gsize (2))] ; 

X coefficient for z*“(l-l) is located at column (n+l-(l-l)), 

X additional shift left because of {}z ! 

end 

GF_num( : , i)=GF_num( : , l)+Gamma , *FProd*Fi( : , ' ,n)*F2 ( : , : ,n) 

* [zeros (Gsize (2) , Gsize (2) ) ; eye (Gsize (2) )] ; 

X since z_(k+n)*0 

GF_num( : ,n)«GF_num(: ,n)+F5* [eye(Gsize (2)> ; zeros (Gsize (2) ,Gsize(2))] ; 

GF_num(: ,n+l)*F5* [zeros (Gsize (2) ,Gsize(2)) ; eye (Gsize (2))] ; 


A. 4 Hinf-calc.m 

X This m-file calculates the H-infinity norms of various transfer functions with the 
X gust measurement data as input "u" and the particular variable as the output H y w . 

X For a discrete-time state space representation I can make use of the Matlab 
X function normhinf . Like for the rest of the simulation, I have to restrict to 
X scalar gust data. The dynamic matrix Aall refers to the controlled system, Atall 
X to the case of no controller in use 

L*RC/norm (RC , r f ro 1 ) ; 

Bi_mat=zeros(Bsize(2) ,n+l) ; 
if exist ( ’Thetat* , , var , )*=l 
for 1*1 :n+l 

Bi_mat ( : ,l)=F3*GF_nuin( : ,n+2-l) ; X for approach II 

end 

else 
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for 1*2 :n+l 

Bi_mat ( : ,l)«F3*GF_num( : ,n+2-l) ; 7* for approach I 

end 

end 

K*F3*F4 ; 

All*Phi+Gamma*K ; 7. Phi_c 

A12“Gamma*Bi_mat+ [zeros (Gsize) Theta zeros (Gsize (1) , (n-l)*Gsize (2) )] ; 

A21*zeros(n+l .Asize (1)) ; 

A22*zeros(n+l ,n+l) ; 
for l-i:n 

A22(i,i+i)-i; 

end 

Aall-CCAli A12] ; [A21 A22]] ; 

At 11 'Phi; 

Atl2“ [zeros (Gsize) Theta zeros(Gsize (1) , (n-l)*Gsize (2) )] ; 

At 21 “A 21 ; 

At22«A22; 

Atall* [[Atll Atl2] ; [At21 At22]]; 

Ball*zeros (Asize (l)+n+l , 1) ; 

Ball (Asize (l)+n+l , l)“i ; *1 last element of column vector ie 1 

Callx* [eye (Asize (1) ) zeros (Asize (1) ,n+l)] ; '/, for the "system" m->system->x 
Callu-[K Bi^mat] ; % for the "system" m->system->u 

if exist (’Thet at’ . , var , )--l % approach II ? 

if axiatC’Ml’ , ’var’)**! '/, PT2 actuator ? 

Calljn*-L*inv(M2)* [K2 D2 zeros (D2size(l) ,2*Dlsize(2) )3 

* [A zerosCGsize) G zeros (Gsize(l) , (n-l)*Gsize (2))] ; 

else 

Calljn— L*inv(M2)* [K2 D2 zeros (D2size(l) ,Dlsize(2) )] 

* [A zeros(Gsize) G zeros (Gsize(l) , (n-l)*Gsize (2) )] ; 

end 

Calljc*L*( [R* (A+B*K) R*B*Bi_mat] 

+ [zeros (D2size(l) ,Asize(2)) Gt* [zeros (Gsize (2) ,Gsize(2)) ;eye(Gsize(2))] 

Gt* [eye (Gsize (2)) ; zerosCGsize (2) , Gsize (2) )] zeros (D2size(l) , (n-l)*Gsize (2) )] ) ; 
else 7* approach I 

Callan*[L*A zeros(Gsize) L*G zeros (Gsize(l) , (n”l)*Gsize(2))] ; 

Callac*L*( [(A+B*K) B*Bi_mat] 

+ [zeros(Asize) zeros(Gsize) G zeros (Gsize(l) , (n-l)*Gsize(2) )]) ; 

end 

Dallx*zeros (Asize (2) , 1) ; 

Dallu-zeros (Bsize (2) , 1) ; 

Dallan*zeros(Asize(2) , 1) ; 

Dallac-Dallan; 

Dalljn*zeros(D2size(2) ,1) ; 

Dalljc-Dalljn; 

sysx=ss (Aall , Ball , Callx , Dallx , T) ; 
sysu*ss (Aall , Ball , Callu , Dallu , T) ; 
if exist (’Thetat* , ’varO"! 

sy8jn*ss(Atall .Ball, Call jn.Dalljn.T) ; 
sys j c»bb (A all , Ball , Call j c , Dali j c , T) ; 
else 

sysan=ss (Atall , Ball , Callan , Dalian , T) ; 
sysac'ss (Aall .Ball, Callac,Dallac,T) ; 

end 


7. approach II ? 
7. approach I 
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nonnx*normhinf (sysx) 
normu*normhinf (sysu) 

if existC’Thetat * , ‘var’)"! */. approach II ? 

normjn*normhinf (sysjn) 
normjc*normhinf (sysjc) 
improve, j- (norm jn-normjc)/normjc 
else % approach I 

norman-normhinf (sysan) 
normac*normhinf (sysac) 
improve (norman-normac ) /normac 

end 


A. 5 pzmap-calc.m 

'/, This m-filo calculates the Pole-zero map of the system with and without 
'/, the controller 

A=Phi ; 

Ac»Phi+Gamma*F3*F4 ; 
if exist ( ’Ml* , ’var’ )“1 

O [eye (D2size) , zeros (D2size) , zeros(D2size(l) ,2*Dlsize(2))] ; 

else 

C*[eye (D2size) , zeros (D2size) , zeros(D2size (1) ,Dlsize(2))] ; 

end 

B-zeros (Bsize) ; */, B, D do not matter for this purpose 

D-zeros CD2size (1) ,Bsize(2)) ; 

sys«BS (A ,B,C, D,T) ; 
sysc^ss (Ac,B, C,D ,T) ; 

[P , Z] ■'pzmap ( sy s ) ; 

[Pc ,Zc] -pzmap(sysc) ; 

dummy*tf (1 , 1) ; 

pzmap (dummy) ; '/, just to get the lables etc. 

zgrid; 

hold on 

plot (P , ’gx’); 

plot (Pc, ’bx’); 

SZ=size (Z) ; 

SZc*=size (Zc) ; 
if SZ(l)~-0 

plot (Z, J go>); 

end 

if SZc(l)--0 

plot(Zc, J bo'); 

end 


